# Manuscript: Signaling Resolve Through Credit-Claiming
# Author: Ilayda B. Onder, The Pennsylvania State University, Department of Political Science
# Replication File for International Interactions
# April 8, 2023

# This analysis is run using R version 4.2.1 (2022-06-23) for Mac
# The codes in this replication file produces Table A2, Table A3, Table A4, Table A5, Table A6, Table A7, Table A8, Table A9, Table A10, Table A11, Table A12, Table A13, Table A14, Table A15, Table A16 presented in the Online Appendices
# The sections in this replication file are in order. Please run them one by one. 

# The packages needed to be installed to run the analysis:
install.packages("stargazer") # for publication quality tables
install.packages("sampleSelection") # for Heckman selection models
install.packages("pscl") # for Zero-inflated negative binomial models
install.packages("systemfit") # for seemingly unrelated regression models
install.packages("ggplot2") # for graphs
install.packages("dplyr") # for dataframe manipulation
install.packages("texreg") # for tables

# Load packages
library(stargazer)
library(sampleSelection)
library(pscl)
library(systemfit)
library(ggplot2)
library(dplyr)
library(texreg)

# WARNING: Please set the working directory to the current Source File Location by clicking Session > Set Working Directory > To Source File Location

# Load analysis datasets
load("Analysis Data Part 2.RData")

dt$num_rivals = as.numeric(as.character(dt$num_rivals)) # convert number of rivals variable to numeric
dt$unclaimed_log = log(dt$unclaimed + 1) # create the logged version of the unclaimed attacks variable
dt$all = dt$claimed + dt$unclaimed # create total number of attacks variable using unclaimed and claimed attacks variables
dt$all_log = log(dt$all + 1) # create the logged version of the total number of attacks variable
dt$claimedratio = dt$claimed / dt$all # create ratio of claimed attacks to the total number of attacks
dt$rivaldum = ifelse(dt$num_rivals == 0, 0, 1) # create dummy variable for the presence of rival groups using the number of rivals
dt$enemydum = ifelse(dt$num_enemies == 0, 0, 1) # create dummy variable for the presence of enemy groups using the number of enemies


# Appendix 3: Table A2: Table 2's zero components----
## Selection models: 1st stage is whether there was at least one attack, 2nd stage is whether all attacks were claimed
dt$stage1 = ifelse(dt$all > 0, 1, 0)
dt$stage2 = ifelse(dt$claimedratio == 1, 1, 0)

summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi,
                        data = dt))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = dt))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = dt))

summary(s4 <- selection(stage1 ~  + 
                          mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = dt))

summary(s5 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi +
                          lead_hierarch + rivaldum,
                        data = dt))

summary(s6 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = dt))

summary(s7 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = dt))

summary(s8 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = dt))

stargazer(s1, s2, s3, s4, s5, s6, s7, s8, selection.equation = TRUE,# only print the selection equation
          type = "text", 
          covariate.labels = c("Multiple base countries",
                               "Group duration (in years)",
                               "Decapitation (l1)",
                               "Constant"),
          dep.var.labels = c("DV: Whether the group perpetrated an attack"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 4: Table A3: OLS: DV is the percentage of claimed attacks----
summary(o0 <- lm(claimed_per ~
                   duration + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi,
                 data = dt))

summary(o1 <- lm(claimed_per ~
                   duration*splinterhist +
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel,
                 data = dt))

summary(o2 <- lm(claimed_per ~
                   duration*civi + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel,
                 data = dt))

summary(o3 <- lm(claimed_per ~
                   fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi,
                 data = dt))

summary(o4 <- lm(claimed_per ~
                   duration + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o5 <- lm(claimed_per ~
                   duration*splinterhist +
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o6 <- lm(claimed_per ~
                   duration*civi + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o7 <- lm(claimed_per ~
                   fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi +
                   lead_hierarch + rivaldum,
                 data = dt))

stargazer(o0, o1, o2, o3, o4, o5, o6, o7, 
          type = "text",
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          dep.var.labels = c("DV: Percentage of Claimed Attacks"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 4: Table A4: OLS: DV is the percentage of claimed attacks (drop zero attack years)----
summary(o0 <- lm(claimedratio ~
                   duration + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi,
                 data = dt))

summary(o1 <- lm(claimedratio ~
                   duration*splinterhist +
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel,
                 data = dt))

summary(o2 <- lm(claimedratio ~
                   duration*civi + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel,
                 data = dt))

summary(o3 <- lm(claimedratio ~
                   fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi,
                 data = dt))

summary(o4 <- lm(claimedratio ~
                   duration + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o5 <- lm(claimedratio ~
                   duration*splinterhist +
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o6 <- lm(claimedratio ~
                   duration*civi + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel +
                   lead_hierarch + rivaldum,
                 data = dt))

summary(o7 <- lm(claimedratio ~
                   fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                   all_log + 
                   nonterr_casualties_log + diversity + shr_trans + rel + 
                   splinterhist + civi +
                   lead_hierarch + rivaldum,
                 data = dt))

stargazer(o0, o1, o2, o3, o4, o5, o6, o7, 
          type = "text",
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          dep.var.labels = c("DV: Percentage of Claimed Attacks"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)



# Appendix 5: Table A5: Seemingly unrelated regression----
summary(su0 <- systemfit(list(claimed = claimed_log ~ duration +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi, 
                              unclaimed = unclaimed_log ~ duration +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi), 
                         data = dt, method = "SUR"))

summary(su1 <- systemfit(list(claimed = claimed_log ~ duration*splinterhist +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel, 
                              unclaimed = unclaimed_log ~ duration*splinterhist +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel), 
                         data = dt, method = "SUR"))

summary(su2 <- systemfit(list(claimed = claimed_log ~ duration*civi +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel, 
                              unclaimed = unclaimed_log ~ duration*civi +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel), 
                         data = dt, method = "SUR"))

summary(su3 <- systemfit(list(claimed = claimed_log ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi, 
                              unclaimed = unclaimed_log ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi), 
                         data = dt, method = "SUR"))

summary(su4 <- systemfit(list(claimed = claimed_log ~ duration +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi +
                                lead_hierarch + rivaldum, 
                              unclaimed = unclaimed_log ~ duration +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi +
                                lead_hierarch + rivaldum), 
                         data = dt, method = "SUR"))

summary(su5 <- systemfit(list(claimed = claimed_log ~ duration*splinterhist +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel +
                                lead_hierarch + rivaldum, 
                              unclaimed = unclaimed_log ~ duration*splinterhist +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel +
                                lead_hierarch + rivaldum), 
                         data = dt, method = "SUR"))

summary(su6 <- systemfit(list(claimed = claimed_log ~ duration*civi +
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel +
                                lead_hierarch + rivaldum, 
                              unclaimed = unclaimed_log ~ duration*civi +
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel +
                                lead_hierarch + rivaldum), 
                         data = dt, method = "SUR"))

summary(su7 <- systemfit(list(claimed = claimed_log ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                                unclaimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi +
                                lead_hierarch + rivaldum, 
                              unclaimed = unclaimed_log ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                                claimed_log +
                                nonterr_casualties_log + diversity + shr_trans + rel + 
                                splinterhist + civi +
                                lead_hierarch + rivaldum), 
                         data = dt, method = "SUR"))

# WARNING: screenreg command produces the table in a messy form. I have manually re-arranged the rows in Word, please refer to the variable names when replicating Table A5
texreg::screenreg(list(su0, su1, su2, su3, su4, su5, su6, su7))


# Appendix 6: Table A6: Selection models with group goal controls----
## Selection models: 1st stage is whether there was at least one attack, 2nd stage is whether all attacks were claimed
dt$stage1 = ifelse(dt$all > 0, 1, 0)
dt$stage2 = ifelse(dt$claimedratio == 1, 1, 0)

summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi +
                          ercsr + tch,
                        data = dt))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          ercsr + tch,
                        data = dt))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          ercsr + tch,
                        data = dt))

summary(s4 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          ercsr + tch,
                        data = dt))


stargazer(s1, s2, s3, s4, type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "ercsr", "tch"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Regime-change goal",
                               "Territorial-change goal",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)


# Appendix 6: Table A7: Zero-inflated negative binomial models with group goal controls----
summary(z1 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi +
                         ercsr + tch |
                         duration,
                       data = dt, dist = "negbin"))

summary(z2 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         ercsr + tch |
                         duration,
                       data = dt, dist = "negbin"))

summary(z3 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         ercsr + tch |
                         duration,
                       data = dt, dist = "negbin"))

summary(z4 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         ercsr + tch |
                         fate_leader_dummy,
                       data = dt, dist = "negbin"))

stargazer(z1, z2, z3, z4, type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "ercsr", "tch"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Regime-change goal",
                               "Territorial-change goal",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 6: Table A8: Selection models with non-violent group activities controls----
dt$social_service[dt$social_service == 9] <- NA # 9 indicates a missing observation, replace it with NA
summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi +
                          state_sponsor + public_service + social_service + drugtk,
                        data = dt))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist+ 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          state_sponsor + public_service + social_service + drugtk,
                        data = dt))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          state_sponsor + public_service + social_service + drugtk,
                        data = dt))

summary(s4 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          state_sponsor + public_service + social_service + drugtk,
                        data = dt))


stargazer(s1, s2, s3, s4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "state_sponsor", "public_service", "social_service", "drugtk"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "State sponsorship",
                               "Public service provision",
                               "Social service provision",
                               "Drug-trafficking",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)


# Appendix 6: Table A9: Zero-inflated negative binomial models with non-violent group activities controls----
summary(z1 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi +
                         state_sponsor + public_service + social_service + drugtk |
                         duration,
                       data = dt, dist = "negbin"))

summary(z2 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         state_sponsor + public_service + social_service + drugtk |
                         duration,
                       data = dt, dist = "negbin"))

summary(z3 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         state_sponsor + public_service + social_service + drugtk |
                         duration,
                       data = dt, dist = "negbin"))

summary(z4 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         state_sponsor + public_service + social_service + drugtk |
                         fate_leader_dummy,
                       data = dt, dist = "negbin"))
stargazer(z1, z2, z3, z4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "state_sponsor", "public_service", "social_service", "drugtk"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "State sponsorship",
                               "Public service provision",
                               "Social service provision",
                               "Drug-trafficking",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 6: Table A10: Selection models with regional controls----
summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi +
                          EAP + ECA + LAC + MENA + SAS + SSA,
                        data = dt))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist+ 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          EAP + ECA + LAC + MENA + SAS + SSA,
                        data = dt))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          EAP + ECA + LAC + MENA + SAS + SSA,
                        data = dt))

summary(s4 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          EAP + ECA + LAC + MENA + SAS + SSA,
                        data = dt))


stargazer(s1, s2, s3, s4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "EAP", "ECA", "LAC", "MENA", "SAS", "SSA"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "East Asia and Pacific",
                               "Europe and Central Asia",
                               "Latin America and Caribbean",
                               "Middle East and North Africa",
                               "South Asia",
                               "Sub-Saharan Africa",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)


# Appendix 6: Table A11: Zero-inflated negative binomial models with regional controls----
summary(z1 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi +
                         EAP + ECA + LAC + MENA + SAS + SSA |
                         duration,
                       data = dt, dist = "negbin"))

summary(z2 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         EAP + ECA + LAC + MENA + SAS + SSA |
                         duration,
                       data = dt, dist = "negbin"))

summary(z3 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         EAP + ECA + LAC + MENA + SAS + SSA |
                         duration,
                       data = dt, dist = "negbin"))

summary(z4 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         EAP + ECA + LAC + MENA + SAS + SSA |
                         fate_leader_dummy,
                       data = dt, dist = "negbin"))
stargazer(z1, z2, z3, z4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "EAP", "ECA", "LAC", "MENA", "SAS", "SSA"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "East Asia and Pacific",
                               "Europe and Central Asia",
                               "Latin America and Caribbean",
                               "Middle East and North Africa",
                               "South Asia",
                               "Sub-Saharan Africa",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 7: Table A12: Only groups that use suicide terrorism: Selection models----
dt$stage1 = ifelse(dt$all > 0, 1, 0)
dt$stage2 = ifelse(dt$claimedratio == 1, 1, 0)

summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi,
                        data = subset(dt, suic == 1)))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, suic == 1)))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, suic == 1)))

summary(s4 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, suic == 1)))

summary(s5 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi +
                          lead_hierarch + rivaldum,
                        data = subset(dt, suic == 1)))

summary(s6 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = subset(dt, suic == 1)))

summary(s7 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = subset(dt, suic == 1)))

summary(s8 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel +
                          lead_hierarch + rivaldum,
                        data = subset(dt, suic == 1)))

stargazer(s1, s2, s3, s4, s5, s6, s7, s8, type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 7: Table A13: Only groups that use suicide terrorism: Zero-inflated negative binomial models----
summary(z1 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z2 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z3 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z4 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         fate_leader_dummy,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z5 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi +
                         lead_hierarch + rivaldum |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z6 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z7 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         duration,
                       data = subset(dt, suic == 1), dist = "negbin"))

summary(z8 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         fate_leader_dummy,
                       data = subset(dt, suic == 1), dist = "negbin"))

stargazer(z1, z2, z3, z4, z5, z6, z7, z8, type = "text",
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 7: Table A14: Only groups with short lifespans: Selection models----

last_row <- dt %>% # select the last row for each group, last row indicates the maximum group duration
  group_by(gname) %>%
  slice(n())
last_row = last_row[c("gname", "duration")] # only keep the group name and duration variables
colnames(last_row)[2] <- "lifespan" # change column names
summary(last_row$lifespan) # maximum group duration is 47 years
dt = merge(dt, last_row, by = "gname", all.x = T) # merge maximum group duration data with the main dataframe
dt$short = ifelse(dt$lifespan > 23.5, 0, 1) # create a dummy variable for groups with short lifespans

dt$stage1 = ifelse(dt$all > 0, 1, 0)
dt$stage2 = ifelse(dt$claimedratio == 1, 1, 0)

summary(s1 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel + 
                          splinterhist + civi,
                        data = subset(dt, short == 1)))

summary(s2 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*splinterhist + 
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, short == 1)))

summary(s3 <- selection(stage1 ~ mul_bases +
                          duration,
                        stage2 ~ 
                          duration*civi +
                          all_log +
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, short == 1)))

summary(s4 <- selection(stage1 ~ mul_bases + 
                          fate_leader_dummy,
                        stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                          splinterhist + civi +
                          all_log + 
                          nonterr_casualties_log + diversity + shr_trans + rel,
                        data = subset(dt, short == 1)))


stargazer(s1, s2, s3, s4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)

# Appendix 7: Table A15: Only groups with short lifespans: Zero-inflated negative binomial models----


summary(z1 <- zeroinfl(claimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi |
                         duration,
                       data = subset(dt, short == 1), dist = "negbin"))

summary(z2 <- zeroinfl(claimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = subset(dt, short == 1), dist = "negbin"))

summary(z3 <- zeroinfl(claimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = subset(dt, short == 1), dist = "negbin"))

summary(z4 <- zeroinfl(claimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         fate_leader_dummy,
                       data = subset(dt, short == 1), dist = "negbin"))

stargazer(z1, z2, z3, z4, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)


# Appendix 8: Table A16: Zero-inflated negative binomial models of unclaimed attacks----
summary(z1 <- zeroinfl(unclaimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi |
                         duration,
                       data = dt, dist = "negbin"))

summary(z2 <- zeroinfl(unclaimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = dt, dist = "negbin"))

summary(z3 <- zeroinfl(unclaimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         duration,
                       data = dt, dist = "negbin"))

summary(z4 <- zeroinfl(unclaimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel |
                         fate_leader_dummy,
                       data = dt, dist = "negbin"))

summary(z5 <- zeroinfl(unclaimed ~ 
                         duration + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel + 
                         splinterhist + civi +
                         lead_hierarch + rivaldum |
                         duration,
                       data = dt, dist = "negbin"))

summary(z6 <- zeroinfl(unclaimed ~ 
                         duration*splinterhist + 
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         duration,
                       data = dt, dist = "negbin"))

summary(z7 <- zeroinfl(unclaimed ~ 
                         duration*civi +
                         all_log +
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         duration,
                       data = dt, dist = "negbin"))

summary(z8 <- zeroinfl(unclaimed ~ 
                         fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                         splinterhist + civi +
                         all_log + 
                         nonterr_casualties_log + diversity + shr_trans + rel +
                         lead_hierarch + rivaldum |
                         fate_leader_dummy,
                       data = dt, dist = "negbin"))
stargazer(z1, z2, z3, z4, z5, z6, z7, z8, type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)",
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"),
          model.numbers = T)
# END----

